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ABSTRACT 

Local Optima Networks (LONs) have been recently proposed as an 
alternative model of combinatorial fitness landscapes. The model 
compresses the information given by the whole search space into 
a smaller mathematical object that is the graph having as vertices 
the local optima and as edges the possible weighted transitions be- 
tween them. A new set of metrics can be derived from this model 
that capture the distribution and connectivity of the local optima in 
the underlying configuration space.This paper departs from the de- 
scriptive analysis of local optima networks, and actively studies the 
correlation between network features and the performance of a lo- 
cal search heuristic. The NK family of landscapes and the Iterated 
Local Search metaheuristic are considered. With a statistically- 
sound approach based on multiple linear regression, it is shown 
that some LONs' features strongly influence and can even partly 
predict the performance of a heuristic search algorithm. This study 
validates the expressive power of LONs as a model of combinato- 
rial fitness landscapes. 

Categories and Subject Descriptors 

F.2.m [Analysis of Algorithms and Problem Complexity]: Mis- 
cellaneous; G.2.2 [Discrete Mathematics]: Graph Theory — Net- 
work problems; 1.2.8 [Artificial Intelligence]: Problem Solving, 
Control Methods, and Search — Heuristic methods 

General Terms 

Algorithms, Measurement, Performance 

Keywords 

Combinatorial Fitness Landscape, Local Optima Network, Local 
Search Heuristics 
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1. INTRODUCTION 

One of the most conspicuous limitations of heuristic search meth- 
ods in combinatorial optimization, is the ability to become trapped 
at a local optimum |9 |. The number and distribution of local op- 
tima in a search space have, therefore, an important impact on the 
performance of heuristic search algorithms. A recently proposed 
model: Local Optima Networks |22| |25| , provides an intermedi- 
ate level of description for combinatorial fitness landscapes. The 
model has a higher descriptive power than a single statistical metric 
or set of metrics; but it also compresses the search space into a more 
manageable mathematical object. Specifically, a graph having as 
vertices the optima configurations of the problem and as edges the 
possible weighted transitions between these optima. This network 
representation allows the application of new analytical tools and 
metrics to study combinatorial landscapes, namely, those of com- 
plex networks analysis (e.g. degree distribution, clustering coeffi- 
cient, assortativity, and community structure, to name a few). In 
previous work, alternative definitions of edges have been studied, 
and some of these metrics have been computed on the network ex- 
tracted for two combinatorial problems: the NK family of land- 
scapes |22| 1 1 6| [T71 125 1 124) , and the Quadratic Assignment prob- 
lem J7]|6j. Those studies have been mainly descriptive, although 
distinctive correlations between some network features and previ- 
ous knowledge about search difficulty in these landscapes had been 
found. A previous related work [18], explored the relationships be- 
tween NK landscape features and the performance of a hybrid EA. 
The authors use standard landscape metrics, and conduct a study 
based mainly on scatter plots. They suggest that: "further work is 
necessary to gain better understanding of the escape rate on the ac- 
tual problem difficulty". The present study, addresses exactly this 
point. The goal is to systematically explore correlations between 
local optima network features and the performance of a stochas- 
tic local search algorithm (Iterated Local Search) running on the 
underlying combinatorial optimization problem (in this study the 
NK family of landscapes). The 'escape rate' is a property re- 
lated to the local optimum, and the LON could be considered as a 
new tool to better understand problem difficulty. The ultimate goal 
is to have predictive models of the performance of specific local 
search heuristics when solving a given combinatorial optimization 
problem. This paper proposes an initial predictive model of perfor- 
mance based on the most influential LON features. 



2. METHODS 

The local optima network model of combinatorial landscapes, 
and the iterated local search metaheuristic are considered in this 
study. The relevant definitions and experimental setup are given 
below. 

2.1 Local Optima Network 

A fitness landscape (20) is a triplet (S, V, f) where S is a set 
of potential solutions i.e. the search space, V : S — > 2 s is a 
function that assigns to every s£ S a set of neighbors V(s) i.e the 
neighborhood structure, and / : S — > M is the evaluation of the 
corresponding solutions i.e. the fitness function. 

The present study uses the well-known A*' if -landscapes |12| as 
a benchmark set. It is a problem-independent model for construct- 
ing combinatorial landscapes that are tunably rugged. In the model, 
N refers to the number of (binary) genes in the genotype, i.e. the 
string length, and K to their epistatic interaction, i.e. the number 
of other loci (chosen at random here) that influence the fitness con- 
tribution of a particular gene. Starting from this iV-loci 2-allele 
additive model, by increasing the non-linearity K from to N — 1, 
the landscapes can be tuned from smooth to rugged. Hence, S is 
the search space of all iV-bit binary strings, and its size is $S = 2 . 
The neighborhood is defined by the minimum possible move on it, 
which is the single bit-flip operation, and the neighborhood size 
is jj V(s) = N. The fitness function evaluates each genotype s 
as f(s) = J2iLi fi( s i> sn,. . . , s iK ), where the values of loci 
contributions /; : {0, — > [0, 1] are drawn uniformly at ran- 
dom in [0, 1]. 

A local optimum (LO), which is taken to be a maximum here, is 
a solution s* such that Vs G V(s), f(s) < f(s*). All optima are 
determined through exhaustive search by recursively running the 
1 -bit-flip best-improvement hill-climber, as in Algorithm[T] Let us 
denote by h(s) the operator that associates to each solution s g S, 
the solution obtained after applying that algorithm until conver- 
gence to a LO. Since S is of finite size and there is no neutrality 
in f(s) values, this produces a partition of the landscape in a finite 
number of basins of attraction, which we can denote by LO\, LOi, 
LO3 . . . , LO„ v , the local optima. 



Algorithm 1: Best-Improvement Hill-Climber 

Choose initial solution s G S ; 
repeat 

choose s' € V(s), such that f(s') = max xeV (s)f(x); 
if/0) </(s')then 

until s is a Local Optimum; 



The connections among them account for the chances of escap- 
ing from a LO and jumping into another one with a controlled 
move [24|. There exists a directed transition ey from LOi to LOj 
if it exists a solution s such that d(s,LOi) < D and h(s) = 
LOj, where the distance d(si,Sj) can be measured in "number 
of moves" (i.e Hamming distance in the bit-flip operator case). The 
distance-threshold D G N can be chosen accordingly to the ap- 
plied perturbation; in this work, it is set to D = 2. The weight 
Wij of such a transition is then: Wij = (t{s G S | d(s,LOi) < 
D and h(s) = LOj}, i.e. the number of paths at distance D start- 
ing at LOi and reaching the basin of LOj. This can be normalized 
by the number of solutions within reach w.r.t. the given distance 
threshold, i.e. ${s e S \ d{s, LOi) < D}. 

The weighted and directed graph G = (V, E) having the set 



of vertices V = {LO\, . . . , LO„ v } and the set of edges E = 
{eij \wij > 0}, is the Local Optima Network (LON) (24) . 

2.2 Iterated Local Search 

Iterated local search is a relatively simple but successful algo- 
rithm. It operates by iteratively alternating between applying a 
move operator to the incumbent solution and restarting local search 
from the perturbed solution. This search principle has been redis- 
covered multiple times, within different research communities and 
with different names |3. 14 1. The term iterated local search (ILS) 
was proposed in |13| . Algorithm[2]outlines the procedure. 



Algorithm 2: Iterated Local Search 

so GeneratelnitialSolution; 
s* LocalSearch(so); 
repeat 

s' <— Perturbation(s*); 

a'* i— LocalSearch(s ); 
s* <s— AcceptanceCriterion(s*, s'*); 
until termination condition met; 



In the present study, the base LocalSearch heuristic is the same 
best-improvement hill-climber of Algorithm [T] which stops on a 
LO. This heuristic uses the single bit-flip move operator. Therefore, 
a 2-bit-flip mutation is chosen as a Perturbation operator. When a 
different LO is found after that, the search process accepts the move 
if its fitness is higher (we are assuming maximization). With these 
settings, ILS is performing a first-improvement hill-climbing in the 
configuration space of the LON with the escape-edges at distance 
D = 2, as defined in section [2~T1 

The search terminates at the global optimum, which for bench- 
mark problems is known a priori, or when reaching a pre-set limit 
of fitness evaluations FE max . 

2.3 Performance Evaluation 

As the performance criterion, we selected the expected number 
of function evaluations to reach the global optimum (success) after 
independent restarts of the ILS algorithm (Algorithm|2| 1 1 1. This 
measure accounts for both the success rate (p s G (0, 1]) and the 
convergence speed. In theory, after (N — 1) unsuccessful runs 
stopped at T us -steps and the final successful one running for T s - 
steps, the total running time would be T = X^^Ji 1 {T us ) k + T s . 
Taking the expectation and considering that TV follows a geometric 
distributiorQwith parameter p s , it gives: 

E(T) - (^y^J E(T U6 ) + E(T S ) 

where in the present case E(T US ) = FE max , the ratio of successful 
to total runs is an estimator for p s , and E(T a ) can be estimated by 
the average running time of those successful runs. 

The ILS variant detailed in Sec. |2.2| is essentially incomplete, 
i.e. there are soluble problem instances for which the success prob- 
ability is < 1 even in the limit of an infinite running time |1 1 1 . 
Given the chosen acceptance criterion, the search will eventually 
get stuck. Indeed, out of the test runs, 1 instance with K = 16 and 
3 instances with K = 17 were not solved. This theoretical limita- 
tion could be overcome by performing as many random restart as 
to cover the whole search space, but such a solution is of limited 

'probability distribution of the number N— 1 of failures before the 
first success (Bernoulli trials) 



practical interest for large problems. In the present study, the suc- 
cess performance has been estimated on all the instances that had 
been solved at least once. 

The benchmark set consists of TVif-landscapes with TV = 18 
and K £ {2, 4, 6, 8, 10, 12, 14, 16, 17}. Those are the largest pos- 
sible parameter combinations for which we could afford the ex- 
haustive extraction of the local optima networks. In order to mini- 
mize the influence of the random generation of landscapes, 30 in- 
dependent problem instances are considered for each combination 
of TV and K, which accounts for a total of 270 instances in the 
problem set. The function-evaluations limit is set to 1/5 of JjS, i.e. 
FEmax — 5.2 • 10 4 , success rate p s and running time of successful 
runs T s are estimated on 500 random restarts per instance. 

3. RESULTS 

3.1 Descriptive Statistics 

Table [T] summarises the LON metrics. Results are grouped ac- 
cording to the K value of the corresponding TV if -landscapes, and 
present averages and standard deviations over 30 independent re- 
alizations per group. The number of local optima nv, is a metric 
familiar to any description of a rugged landscape. The other met- 
rics are particular to the complex-network perspective provided by 
the LON model Q5). 

From left to right in the table, lo represents the average length 
of the shortest paths that reach the global optimum starting from 
any other local optimum. The cost associated to an edge ij of the 
LON graph, is dkj = l/wtj. This measure can be interpreted as 
the expected number of random perturbations for escaping LOi 
and entering exactly the basin of LOj . Iv gives the average path 
length for the whole graph, which accounts for the weighted net- 
work characteristic length, lo is intuitively more directly related 
to the search difficulty. Indeed, lo increases steadily with the land- 
scape ruggedness K, whilst the trend with Iv is less clear. As a 
possible explanation, the network growth in terms of nodes might 
be counteracted by a growth in nodes connectivity: the number of 
weighted outgoing transitions from a given LO, i.e. its out-degree 
in the LON, increases with K (cf. column zout in Tab.[T](. 

Column fnn measures the correlation between the fitness of a 
node and the weighted average of the fitness of its nearest neigh- 
bors (2). This is relevant as the ILS acceptance criterion takes the 
fitness values of LO into account. With respect to this metric, NK- 
landscapes behave at the LO level as they do at the solution level: 
they become the more and more uncorrected with K approaching 
TV — 1 (see column fnn). In general, it is expected that a high and 
positive fnn correlation would help the search process. 

Column wii reports the average number of perturbations that re- 
main in the same basin of attraction, which is a proxy for the basin 
size. The larger the basins, the more difficult it is to escape them. 
Table [T] shows that average wii value decreases with the landscape 
ruggedness. 

Column cc, reports the clustering coefficient |15| , which mea- 
sures the ratio of connected triples in the LON graph. In a social 
network, this coefficient measures the probability for one's friends 
to be friends among each other. In the LO networks, it provides 
an index of topological locality for the transitions between local 
optima. Table [T] suggest that cc values decrease steadily with in- 
creasing K, but so does the LON density |24 | and that might be the 
reason. 

Column zout counts the number of transitions departing from 
a given LO. It is relevant to know whether all the transitions have 
the same rate, or if there is a preferred direction. To this aim, the 
disparity score y2 gauges the weight heterogeneity of out-going 



edges [2|. When all connections ey leaving a given LOi have the 
same probability Wij, the disparity will be close to the inverse of 
the out-degree 1 / zouti otherwise it will be higher than this value. 
Column j/2 shows that the disparity monotonically decreases with 
increasing K. From the point of view of a metaheuristic dynamic, 
a low disparity means that transitions are almost equiprobable. The 
LON topology does not preferentially guide the search trajectory, 
which makes the search process harder. 

Column knn, reports the nearest-neighbors degree correlation or 
assortativity, a classical description of the mixing pattern of nodes 
in a network. The assortativity measures the affinity to connect 
with high or low-degree neighbors. The LONs of TVif-landscapes 
are strongly disassortative, i.e. LO with few connections tend to 
link to others with many, and conversely. The implications of this 
observation on the search difficulty are worth further investigation. 

The last column in Table[T]gives the values of the success perfor- 
mance indicator E(T) detailed in Sec 2.3 which we abbreviate to 
els in the following. Clearly, the expected running time increases 
with the landscape ruggedness (problem non-linearity) K. In order 
to analyze the inter-correlations between all the studied metrics, 
Figure [T] displays a correlogram of the whole data set. The figure 
depicts any possible pairing among the observed variables with a 
scatter plot in the panel below the diagonal, and the corresponding 
correlation coefficient in the upper panel. Table [T] corresponds to 
the first column of the correlation matrix, i.e. the one showing how 
the different network metrics shape against the epistasis K. 

We are mostly interested in the correlations with the performance 
metric ets. Thus, the most relevant scatter plots are those in the last 
row of Figure[T]together with the respective Pearson coefficients in 
the last column. By inspecting these values, the highest positive 
correlation is with the average length of the shortest path climbing 
to the global optimum (r(ets, lo) ~ 0.52). and with the total num- 
ber of LO (r (ets, nv) ~ 0.50). This observation seems reasonable: 
the more rugged the landscape (i.e. the larger its LON), the higher 
the number of hops through the LON to reach the global optimum, 
and thus the longer the expected running time of a restarting lo- 
cal search. Conversely, the larger the LO basins, and the higher the 
nearest-neighbors fitness-fitness correlation, the shorter the running 
time ((r(ets, wii) ~ r(ets, fnn) ~ —0.40)). The catch here, is 
that those metrics are in turn correlated among themselves (e.g. 
r(fnn, wii) ~ 0.92), which prevents one from drawing relation- 
ships of causality from simple pair-wise correlations. 

The last column plots also suggest a strong non-linearity of the 
performance estimator w.r.t. the LON metrics. In order to further 
analyze these relationships, Figure [2] zooms on the relevant scatter 
plots, and displays the logarithm of ets as a function of the con- 
sidered landscape measures. The log-transformation allows to ap- 
proach linearity, highlighting and confirming the results in the last 
column of 11. Namely, the relationships between ets and all the 
metrics but Iv appear clearly. The picture suggests a positive ex- 
ponential trend with lo and zout, and a negative exponential trend 
with fnn, wii, cc, and y2. In the case of nv, the relation could 
also be close to power-law. 

Since data are far from being normal bivariate, a more robust 
measure of association would be the rank-based Spearman's p statis- 
tic, reported in Table [2] These results complement the visual in- 
spection of scatter plots and confirm the previous observations. 

3.2 Statistical Modeling 

1111, along with 11, address already some of the research ques- 
tions asked in Section [TJ but do not provide an explanatory model 
for the algorithm performance as a function of the landscape fea- 
tures. To this end, we perform a multiple linear regression on the 
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Figure 1: Correlation matrix for all pairs of the observed variables, to be read from the diagonal. The lower panel displays scatter 
plots and smoothing splines for every possible pairing; the upper panel gives the corresponding Pearson correlation coefficient, 
with text size proportional to its absolute value. In the lower panel, the smoothing is done through locally-weighted polynomial 
regressions |5|. In the upper panel, the correlation coefficient is tested against the null hypothesis and the resulting p-value is 
symbolically encoded at the levels of 0.1 ('), 0.05 (*), 0.01 (**), and 0.001 (***). 




Figure 2: Scatter plots of the logarithm of the estimated time to succeed against all measured variables (see labels on x-axis). The 
epistasis value K is treated as a category, thus results are in a box-and-whisker diagram per group (fist plot on the top-left corner). 



Table 1: Group averages of all the observed variables, aggregated by the epistasis K of the corresponding TVi^-landscape. Standard 
deviations are given in subscripts, nv = number of vertices (Local Optima), lo = average shortest path to reach the global optimum, lv 
= average path length (dij = X/wij), fnn = Spearman coefficient for the nearest-neighbors fitness-fitness correlation, wii = average 
non-normalized weight of self -loops, cc global clustering coefficient, zout = average out-going degree, y2 = average weight disparity 
for out-going edges, knn = degree assortativity, ets = estimated time to succeed. 
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Table 2: Spearman's p statistic for the correlation between ets 
and the LON metrics (p-value < 2.2e — 16 for all pairings). 
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data, which has the general form: 

Hi = A> + fta;i,i + /3 2 ^, 2 H h PpXi, p + e t (1) 

where the response variable y in our case would be ets and p dif- 
ferent predictors Xj are to be chosen among the LON metrics; e is 
the usual random noise term. 

The least square regression produces estimates f3j for the fij 
model coefficients; the difference between the predicted values and 
the actual, observed values, are the regression residuals: = 

Vi — iji — Vi — Pixn — /3 2 £i2 — PpXi P . 

The difficulty of this analysis is that there are several possible 
explanatory variables, which are in turn intercorrelated. In con- 
sequence, some of them could have a confounding effect on the 
regression. In general, when confounders are known, measurable, 
and measured, it is a good practice to include them in the model. 
We, therefore, start by fitting the following formula: 

log(eis) = I3 + P!k + /3 2 log(nv) + /3 2 lo^ \-/3 10 knn + e (2) 

where, w.r.t the general expression ??, the response and one of the 
predictors have been log-transformed in order to better approach 
linearity, as seen in ??. Moreover, the variable k is qualitative and 
enters the model as a fixed effect, which translates in one appropri- 
ate dummy variable for each class K = 2, K = 4, . . . , K = 17. 

The summary statistics for this model are reported in Table [3] 
In the table caption, the multiple R 2 represents the proportion of 
variance explained by the linear regression. R 2 would be equal to 
1 if all observed data points were lying on the regression plane. 
When comparing models with a different number of predictors, the 
adjusted R 2 should be used instead. The F statistic is the ratio 
of the variance explained by the parameters in the model, to the 
residual or unexplained variance. The p-value is the probability of 
achieving an F that large under the null hypothesis of no effect (4). 

The estimated coefficient /3j and their estimated standard error 
dj are given in the 2 nd and 3 rd columns, respectively. Their ratio 
is the f-statistic (4 th column) that is used to calculate a p-value for 
the significance of the estimation (last column). 



Table 3: Summary statistics of the linear regression model 
on all variables. Residual standard error: 0.8702 on 248 de- 
grees of freedom (4 observations deleted due to missingness). 
Multiple R-squared: 0.8585, Adjusted R-squared: 0.8488. F- 
statistic: 88.52 on 17 and 248 DF, p-value: < 2.2e - 16. 
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0.49062 


-1.8132 


7.10e 


-02 



In this initial model, the average length of paths to the global 
optimum lo is the only predictor with a regression coefficient that 
is statistically-significant at the 0.05 threshold (/3 Jo = 0.04882, 
p-value = 2.37e - 19). 

Therefore, we proceed to perform a step-wise model selection by 
backward elimination |23|. From the initial formula, at each step 
we compute what change in the fit could be produced by dropping 
each predictor in turn, and then we eliminate the one that minimizes 
the AIC score of the resulting model [21J. By iterating until all 
predictors become significant, we obtain the final model: 

log(ets) = /3q + Pi lo + /3 zou tzout + P V 2y2 + fiknnknn + e (3) 

which is detailed in Table [4] 

Table 4: Summary statistics of the final linear regression model. 
Residual standard error: 0.8751 on 261 degrees of freedom (4 
observations deleted due to missingness). Multiple R-squared: 
0.8494, Adjusted R-squared: 0.8471. F-statistic: 368.1 on 4 
and 261 DF, p-value: < 2.2e - 16. 



summary 


Estimate 


Std. Error 


t value 


Pr(>|t|) 


(Intercept) 


10.3838 


0.58512 


17.75 


9.24e - 47 


lo 


0.0439 


0.00434 


10.11 


1.67e - 20 


zout 


-0.0306 


0.00831 


-3.68 


2.81e - 04 


y2 


-7.2831 


1.63038 


-4.47 


1.18c - 05 


knn 


-0.7457 


0.40501 


-1.84 


6.67c - 02 



This final model is able to explain 84.94% of the variance ob- 
served in log(ets) with a linear regression on four variables that 
are all LON network metrics. Among these metrics, the length of 
the paths to the global optimum, and the weight disparity have the 
highest relative importance 1 10 1. 

Without a check on the model assumptions, this would remain 
an observational study and could not be used to make predictions. 
To this end, a combination of parametric tests (not reported for 
space reasons) provided a positive confirmation [19]. However, 
a visual diagnostics is more informative. In particular, Figure [3] 
helps to assess if the regression residuals follow a normal distribu- 
tion with zero mean and homogeneous variance, whereas Figure [4] 
displays the contributions to the model of each predictor in turn, 
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Figure 4: Component+residual plots for the linear regression 
model. Circle dots are, for each observation, the correspond- 
ing residual error from the regression plus the value fitted by 
one explanatory variable alone, plotted against that variable 
(see labels on x-axis). No significant deviation from linearity 
appears (smooth green line against the red dotted line of the 
partial regression). 



Figure 3: Top: residual plot, to asses the hypothesis of 
zero-mean and constant-variance of the regression residu- 
als (circle dots) around the fitted values (dotted line); no 
visually-significant deviation appears (red smooth line). Bot- 
tom: quantile-quantile comparison of the studentized regres- 
sion residuals (circle dots) against the theoretical quantiles (red 



highlighting possible violations of the linearity hypothesis |8|. All 
assumptions seem acceptable. Therefore, formula ?? could be used 
to make inferences. In other words, formula ?? coefficients can 
be interpreted as conditional expectations for the average change 
in the response when one predictor undergoes a unitary change, 
and all the others remain fix. Since the dependent variable is log- 
transformed, this effect would be multiplicative. 

The only important limitation of the proposed model is multi- 
collinearity: due to the complex intercorrelations among LON met- 
rics (cf. ??), predictors are not really independent. This does not 
invalidate the multiple linear regression analysis, but it inflates the 
variance of its coefficients and makes it harder to disentangle their 
respective contributions [ 8 1. 

4. CONCLUSIONS 

This article explored correlations between local optima network 
features and the performance of a stochastic local search algorithm 
running on the underlying combinatorial optimization problem. The 
NK family of landscapes and the iterated local search meta-heuristic 
were considered. It has been shown in previous work e.g. |22||25| 
that some features of the LON networks are related to theTv A- 
landscapes ruggedness, and thus to problem difficulty. However, no 
statistically testable model was presented. The contribution of this 
study was to investigate, with a statistically-sound approach, which 
features of the LONs have a strong influence on the search perfor- 
mance, expressed as expected running times to success. The results 
obtained through the use of a multiple linear regression model show 
that some LON metrics are more important than others. These are: 
the average length of the shortest paths to the optimum, the average 
out-degree, the average disparity, and the degree assortativity. This 
study confirms and provides significant evidence that LON model- 
ing is a compressed-but-relevant view of the fitness landscape, and 
can be used to understand and predict search difficulty. 

It is worth noticing that some network metrics can be estimated 
without knowing the global optimum beforehand, such as the aver- 
age out-degree, the fitness-fitness correlation, the average disparity 
and assortativity. Using these metrics and an adequate statistical 
model, as we have done in this work, opens up exciting possibili- 
ties. With standard sampling methods, larger search spaces could 
be studied. Thereafter, using the performance model based on the 
estimated LON metrics, the search heuristic parameters or its op- 
erators can be off-line tuned, or even on-line controlled. We plan 
to continue working in this direction and to extend this analysis to 
other combinatorial problems such as QAP. 
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